home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / rayview / spheregen.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  6.9 KB  |  303 lines
  1. /*
  2.  * Routines to define a unit sphere Object.
  3.  * If you have a sphere library (-lsphere on most SGI machines), you can
  4.  * use that by defining SPHERELIB on the cc command line.  Otherwise,
  5.  * we use Jon Leech's code to generate a sphere.
  6.  * C. Kolb 6/91
  7.  */
  8. #ifdef SPHERELIB
  9. #include <gl.h>
  10.  
  11. /*
  12.  * Define unit sphere Object.  Here, we make use of sgi's spherelib.
  13.  * If not available, undefine SPHERELIB, and the code below will be used.
  14.  */
  15. Object
  16. GLSphereObjectDefine()
  17. {
  18.     int sphere;
  19.  
  20.     sphere = genobj();
  21.     sphobj(sphere);
  22.     return sphere;
  23. }
  24. #else
  25.  
  26. /*
  27.  * The following code has been slightly modified.  The original sphere
  28.  * code is available via anonymous
  29.  * ftp from weedeater.math.yale.edu (130.132.23.17) in pub/sphere.c
  30.  *
  31.  * Modifications include changing print_triangle to draw a GL object,
  32.  * the inclusion of gl.h,
  33.  * the renaming of main() to GLSphereObjectDefine(), the nuking of anything
  34.  * phigs-related, and setting the ccw flag to be true.
  35.  * 
  36.  * C. Kolb 6/91
  37.  */
  38.  
  39. /*% cc -g sphere.c -o sphere -lm
  40.  *
  41.  * sphere - generate a triangle mesh approximating a sphere by
  42.  *  recursive subdivision. First approximation is an octahedron;
  43.  *  each level of refinement increases the number of triangles by
  44.  *  a factor of 4.
  45.  * Level 3 (128 triangles) is a good tradeoff if gouraud
  46.  *  shading is used to render the database.
  47.  *
  48.  * Usage: sphere [level] [-p] [-c]
  49.  *    level is an integer >= 1 setting the recursion level (default 1).
  50.  *    -p causes generation of a PPHIGS format ASCII archive
  51.  *        instead of the default generic output format.
  52.  *    -c causes triangles to be generated with vertices in counterclockwise
  53.  *        order as viewed from the outside in a RHS coordinate system.
  54.  *        The default is clockwise order.
  55.  *
  56.  *  The subroutines print_object() and print_triangle() should
  57.  *  be changed to generate whatever the desired database format is.
  58.  *
  59.  * Jon Leech (leech@cs.unc.edu) 3/24/89
  60.  */
  61. #include <stdio.h>
  62. #include <math.h>
  63.  
  64. #include <gl.h>    
  65.  
  66. typedef struct {
  67.     double  x, y, z;
  68. } point;
  69.  
  70. typedef struct {
  71.     point     pt[3];    /* Vertices of triangle */
  72.     double    area;    /* Unused; might be used for adaptive subdivision */
  73. } triangle;
  74.  
  75. typedef struct {
  76.     int       npoly;    /* # of triangles in object */
  77.     triangle *poly;    /* Triangles */
  78. } object;
  79.  
  80. /* Six equidistant points lying on the unit sphere */
  81. #define XPLUS {  1,  0,  0 }    /*  X */
  82. #define XMIN  { -1,  0,  0 }    /* -X */
  83. #define YPLUS {  0,  1,  0 }    /*  Y */
  84. #define YMIN  {  0, -1,  0 }    /* -Y */
  85. #define ZPLUS {  0,  0,  1 }    /*  Z */
  86. #define ZMIN  {  0,  0, -1 }    /* -Z */
  87.  
  88. /* Vertices of a unit octahedron */
  89. triangle octahedron[] = {
  90.     { { XPLUS, ZPLUS, YPLUS }, 0.0 },
  91.     { { YPLUS, ZPLUS, XMIN  }, 0.0 },
  92.     { { XMIN , ZPLUS, YMIN  }, 0.0 },
  93.     { { YMIN , ZPLUS, XPLUS }, 0.0 },
  94.     { { XPLUS, YPLUS, ZMIN  }, 0.0 },
  95.     { { YPLUS, XMIN , ZMIN  }, 0.0 },
  96.     { { XMIN , YMIN , ZMIN  }, 0.0 },
  97.     { { YMIN , XPLUS, ZMIN  }, 0.0 }
  98. };
  99.  
  100. /* A unit octahedron */
  101. object oct = {
  102.     sizeof(octahedron) / sizeof(octahedron[0]),
  103.     &octahedron[0]
  104. };
  105.  
  106. /* Forward declarations */
  107. point *normalize(/* point *p */);
  108. point *midpoint(/* point *a, point *b */);
  109. void print_object(/* object *obj, int level */);
  110. void print_triangle(/* triangle *t */);
  111.  
  112. /*extern char *malloc(/* unsigned );*/
  113.  
  114. Object
  115. GLSphereObjectDefine(maxlevel)
  116. int maxlevel;
  117. {
  118.     object *old,
  119.        *new;
  120.     int     ccwflag = 1,    /* Reverse vertex order if true */
  121.         i, level;
  122.     Object sph;
  123.  
  124.     makeobj(sph = genobj());
  125.  
  126.     if (ccwflag) {
  127.     /* Reverse order of points in each triangle */
  128.     for (i = 0; i < oct.npoly; i++) {
  129.         point tmp;
  130.               tmp = oct.poly[i].pt[0];
  131.         oct.poly[i].pt[0] = oct.poly[i].pt[2];
  132.         oct.poly[i].pt[2] = tmp;
  133.     }
  134.     }
  135.  
  136.     old = &oct;
  137.  
  138.     /* Subdivide each starting triangle (maxlevel - 1) times */
  139.     for (level = 1; level < maxlevel; level++) {
  140.     /* Allocate a new object */
  141.     new = (object *)malloc(sizeof(object));
  142.     if (new == NULL) {
  143.         fprintf(stderr, "GLSphereObjectDefine: out of memory, level %d\n",
  144.         level);
  145.         exit(1);
  146.     }
  147.     new->npoly = old->npoly * 4;
  148.  
  149.     /* Allocate 4* the number of points in the current approximation */
  150.     new->poly  = (triangle *)malloc(new->npoly * sizeof(triangle));
  151.     if (new->poly == NULL) {
  152.         fprintf(stderr, "GLSphereObjectDefine: out of memory, level %d\n",
  153.         level);
  154.         exit(1);
  155.     }
  156.  
  157.     /* Subdivide each triangle in the old approximation and normalize
  158.      *  the new points thus generated to lie on the surface of the unit
  159.      *  sphere.
  160.      * Each input triangle with vertices labelled [0,1,2] as shown
  161.      *  below will be turned into four new triangles:
  162.      *
  163.      *            Make new points
  164.      *                a = (0+2)/2
  165.      *                b = (0+1)/2
  166.      *                c = (1+2)/2
  167.      *      1
  168.      *     /\        Normalize a, b, c
  169.      *    /  \
  170.      *    b/____\ c        Construct new triangles
  171.      *    /\    /\            [0,b,a]
  172.      *   /    \  /  \            [b,1,c]
  173.      *  /____\/____\        [a,b,c]
  174.      * 0      a    2        [a,c,2]
  175.      */
  176.     for (i = 0; i < old->npoly; i++) {
  177.         triangle
  178.          *oldt = &old->poly[i],
  179.          *newt = &new->poly[i*4];
  180.         point a, b, c;
  181.  
  182.         a = *normalize(midpoint(&oldt->pt[0], &oldt->pt[2]));
  183.         b = *normalize(midpoint(&oldt->pt[0], &oldt->pt[1]));
  184.         c = *normalize(midpoint(&oldt->pt[1], &oldt->pt[2]));
  185.  
  186.         newt->pt[0] = oldt->pt[0];
  187.         newt->pt[1] = b;
  188.         newt->pt[2] = a;
  189.         newt++;
  190.  
  191.         newt->pt[0] = b;
  192.         newt->pt[1] = oldt->pt[1];
  193.         newt->pt[2] = c;
  194.         newt++;
  195.  
  196.         newt->pt[0] = a;
  197.         newt->pt[1] = b;
  198.         newt->pt[2] = c;
  199.         newt++;
  200.  
  201.         newt->pt[0] = a;
  202.         newt->pt[1] = c;
  203.         newt->pt[2] = oldt->pt[2];
  204.     }
  205.  
  206.     if (level > 1) {
  207.         free(old->poly);
  208.         free(old);
  209.     }
  210.  
  211.     /* Continue subdividing new triangles */
  212.     old = new;
  213.     }
  214.  
  215.     /* Print out resulting approximation */
  216.    
  217.     print_object(old, maxlevel);
  218.     closeobj();
  219.     return sph;
  220. }
  221.  
  222. /* Normalize a point p */
  223. point *normalize(p)
  224. point *p;
  225. {
  226.     static point r;
  227.     double mag;
  228.  
  229.     r = *p;
  230.     mag = r.x * r.x + r.y * r.y + r.z * r.z;
  231.     if (mag != 0.0) {
  232.     mag = 1.0 / sqrt(mag);
  233.     r.x *= mag;
  234.     r.y *= mag;
  235.     r.z *= mag;
  236.     }
  237.  
  238.     return &r;
  239. }
  240.  
  241. /* Return the midpoint on the line between two points */
  242. point *midpoint(a, b)
  243. point *a, *b;
  244. {
  245.     static point r;
  246.  
  247.     r.x = (a->x + b->x) * 0.5;
  248.     r.y = (a->y + b->y) * 0.5;
  249.     r.z = (a->z + b->z) * 0.5;
  250.  
  251.     return &r;
  252. }
  253.  
  254. /* Write out all triangles in an object */
  255. void print_object(obj, level)
  256. object *obj;
  257. int level;
  258. {
  259.     int i;
  260.  
  261.     /* Spit out coordinates for each triangle */
  262.     for (i = 0; i < obj->npoly; i++)
  263.     print_triangle(&obj->poly[i]);
  264. }
  265.  
  266. /* Output a triangle */
  267. void print_triangle(t)
  268. triangle *t;
  269. {
  270.     int i;
  271.     float p[3];
  272.  
  273. #ifdef sgi
  274.     bgnpolygon();
  275. #endif
  276.  
  277.     p[0] = t->pt[0].x; p[1] = t->pt[0].y; p[2] = t->pt[0].z;
  278.     n3f(p);
  279. #ifdef sgi
  280.     v3f(p);
  281. #else
  282.     pmv(p[0], p[1], p[2]);
  283. #endif
  284.     p[0] = t->pt[1].x; p[1] = t->pt[1].y; p[2] = t->pt[1].z;
  285.     n3f(p);
  286. #ifdef sgi
  287.     v3f(p);
  288. #else
  289.     pdr(p[0], p[1], p[2]);
  290. #endif
  291.     p[0] = t->pt[2].x; p[1] = t->pt[2].y; p[2] = t->pt[2].z;
  292.     n3f(p);
  293. #ifdef sgi
  294.     v3f(p);
  295.     endpolygon();
  296. #else
  297.     pdr(p[0], p[1], p[2]);
  298.     pclos();
  299. #endif
  300. }
  301.  
  302. #endif
  303.